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I.  Introduction 


In  radar,  sonar  and  seismology  one  is  frequently  interested  in  estimating  the 
dir ections-of-  arrival  and  the  spectral  densities  of  radiating  sources  from  measure¬ 
ments  provided  by  a  passive  array  of  sensors.  In  the  simplest  case,  the  signals 
received  by  the  sensors  consist  of  scaled  and  delayed  replicas  of  the  waveform 
radiated  by  a  single  source.  In  a  more  complicated  scenario,  there  may  be  multi¬ 
ple  sources  and  multiple  propagation  paths  from  the  sources  to  the  sensors. 

The  problem  of  the  simultaneous  estimation  of  the  directions-of-arrival  and 
the  spectral  densities  of  the  impinging  sources  can  be  regarded  as  a  2-D  (two- 
dimensional)  spectral  estimation  problem.  Given  spatial  and  temporal  samples  of 
the  received  signals,  the  problem  is  to  estimate  the  2-D  spectrum,  or  the  energy 
distribution,  in  both  the  spatial  (wavenumber)  and  temporal  (frequency)  domains, 
i.e,  in  the  wavenumber-frequency  plane.  The  spatial  spectrum  consists  of  point 
masses  at  the  different  wavenumbers.  The  temporal  spectrum  may  consist  of 
point  masses  at  different  frequencies,  in  the  case  of  narrowband  sources,  or  may 
be  continuous  in  the  case  of  wideband  sources.  When  the  spectral-densities  of 
the  sources  are  known,  the  problem  of  direction-of-arrival  estimation  reduces  to  a 
1-D  spectral  estimation  problem  in  the  spatial  domain. 

Since  the  number  of  samples  (i.e.  sensors)  in  the  spatial  domain  is  usually 
small,  classical  Fourier  analysis  yields  low  spatial  resolution.  As  a  result,  alterna¬ 
tive  methods  that  provide  higher  resolution  have  been  developed.  The  Maximum 
Entropy  (ME)  method  of  Burg  (1067)  and  the  Minimum  Variance  (MV)  method 
of  Capon  (1060)  were  the  first  methods  that  provided  increased  resolution  in  the 
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2-D  problem,  while  the  Maximum  Entropy  method  was  originally  developed  for 
the  1-D  problem.  Different  extensions  of  the  Maximum  Entropy  method  to  the 
2-D  problem  have  been  proposed  by  Roucos  and  Childers  (1980),  Lim  and  Malik 
(1981)  and  McClellan  (1982)  among  others.  A  different  high  resolution  method 
for  the  2-D  problem,  known  as  2-D  linear-prediction,  was  developed  by  Jackson 
and  Chicn  (1978),  Jain  and  Ranganath  (1978),  Frost  and  Sullivan  (1979)  and 
Kumareseran  and  Tufts  (1981).  Recently,  another  method  for  the  2-D  problem 
based  on  ARMA  modeling  of  the  received  signals  has  been  proposed  by  Morf  et  al 
(1979),  and  further  developed  and  elaborated  by  Porat  and  Friedlander  (1983) 
and  Nehorai  and  Morf  (1983).  We  should  note  that  the  Maximum-Likelihood 
estimator,  though  asymptotically  optimal,  has  not  been  advocated  for  the  2-D 
problem  because  of  the  high  computational  complexity  involved.  In  the  somewhat 
simpler  1-D  problem,  where  the  spectral  densities  of  the  sources  are  known,  the 
Maximum-Likelihood  estimator  has  been  studied  by  Good  (1963),  Hahn  and 
Tretter  (1973),  Schweppe  (1968)  and  Wax  and  Kailath  (1983a). 

In  the  special  case  that  the  sources  are  narrowband  and  have  the  same 
known  center-frequency,  the  problem  of  direction-of-arrivai  estimation  reduces  to 
the  problem  of  1-D  harmonic-retrieval.  A  method,  tailored  especially  for  this 
problem,  that  provides  better  resolution  than  offered  by  the  Minimum  Variance, 
Maximum  Entropy  and  the  linear- prediction  methods,  was  pioneered  by 
Pisarenko  (1973).  The  method  is  based  on  the  eigenstructure  of  the  covariance 
matrix  of  the  received  signals.  Schmidt  (1979)  (1981),  and  independently  Bien- 
venu  and  Kopp  (1980,  1981),  have  improved  the  resolution  of  Pisarenko's 


method,  and  extended  it  from  the  restricted  uniform  linear  array  case  to  a  gen* 
eral  array.  Schmidt's  method  •  as  Pisarenko’s  -  is  a  time-domain  method  based 
on  the  eigenstructure  of  the  covariance  matrix  of  the  received  signals;  Bienvenu 
and  Kopp’s  method  is  a  frequency-domain  method  based  on  the  eigenstructure  of 
the  spectral  density  matrix  of  the  received  signals.  Related  but  somewhat 
different  eigenstructure  methods  for  the  1-D  problem  have  been  proposed  by  Lig- 
get  (1972),  Owsley  (1973),  Reddi  (1979),  Johnson  and  Degraff  (1982), 
Kumareseran  and  Tufts  (1983)  and  Bronez  and  Cadzow  (1983).  The  methods 
described  above  are  off-line;  on-line  implementations  of  Pisarenko’s  method  have 
been  presented  by  Thompson  (1979),  Cantoni  and  Godara  (1980)  and  Reddy  et 
al.  (1982a).  Extension  of  these  on-line  methods  to  the  2-D  harmonic  retrieval 
problem  have  been  presented  by  Larimore  (1981)  and  Reddy  et.  al  (1982b). 

In  this  paper  we  extend  the  eigenstructure  approach  to  the  2-D  problem. 
We  present  a  time-domain  eigenstructure  method  for  the  2-D  harmonic  retrieval 
problem,  that  is  for  the  simultaneous  estimation  of  both  the  direction-of-arrival 
and  the  center-frequency  of  narrowband  sources,  and  a  frequency-domain  eigen¬ 
structure  method  for  the  estimation  of  the  direction-of-arrival  and  the  spectral 
density  matrix  of  wideband  sources. 

The  narrowband  and  the  wideband  2-D  problems  are  formulated  in  section 
II.  The  time-domain  eigenstructure  method  for  the  narrowband  problem  is 
presented  in  section  ID.  The  frequency-domain  eigenstructure  method  for  the 
wideband  problem  is  presented  in  section  IV.  Simulation  results  that  illustrate 
the  performance  of  the  proposed  algorithms  are  presented  in  section  V. 


II.  Problem  Formulation 


The  formulation  will  be  somewhat  different  for  the  narrowband  and  wide¬ 
band  problems.  Following  the  convention  in  array  processing,  a  problem  will  be 
referred  to  as  a  narrowband  problem  if  the  bandwidth  of  the  impinging  sources  is 
much  smaller  than  the  reciprocal  of  the  propagation  time  of  the  signal  across  the 
array;  otherwise  it  will  be  referred  to  as  wideband. 

A.  The  Narrowband  Problem 

Consider  an  array  with  m  sensors,  each  followed  by  a  tapped-delay-line 
with  p  taps  spaced  D  delay  units  apart.  Assume  that  d  narrowband  sources 
(d  <  mp)  centered  at  frequencies  Uj,...,^,  impinge  on  the  array  from  direc¬ 
tions  0t,...,6d  (see  Fig.  1).  Assume  that  the  signals  emitted  by  the  sources  are 
stationary  zero- mean  narrowband  stochastic  processes.  The  signals  received  by 
the  sensors  are  scaled  and  delayed  replicas  of  the  radiated  signals.  Then,  using 
complex  (analytic)  signal  representation,  the  output  of  the  h- th  delay  unit  of 
the  i-th  sensor,  can  be  expressed  as 

r,(/  -  hD)  =  £  aik8k{t  -  rik  -  hD)  +  n,(<  -  hD)  (1) 

*— i 

where  «*(•)  is  the  signal  radiated  by  the  it-th  source  as  observed  at  an  arbitrarily 
chosen  reference  point,  uk  is  the  center  frequency  of  the  it-th  source,  rik  is  the 
propagation  delay  between  the  «-th  sensor  and  the  reference  point  for  the  it-th 
source,  aik  is  the  amplitude  response  of  the  i-th  sensor  to  the  it-th  source,  and 
n,(  )  is  the  additive  noise  at  the  i-th  sensor. 


Since  «*(•)  is  &  narrowband  process  we  can  well  approximate  the  time-delay 
by  a  phase-shift.  Therefore,  we  can  rewrite  (1)  as 

r,.(!  -  AD)  =  £  +  n,(<  -  AD)  <2> 

Assume  that  the  received  signals  are  sampled  simultaneously  at  times 
tk  (  k  —  ),  yielding  N  "snapshots”,  each  consisting  of  mp  samples 

r,- ( tk  -  hD)  (  i  =  l,...,m;  h  =  0 ,...p  -  1  ).  Grouping  the  samples  correspond¬ 
ing  to  the  i-th  sensor  into  a  p  X  1  vector,  we  can  rewrite  equation  (2)  in  matrix 
form 

«*,(0  =  A  ,-•(/)  +  »,-{*)  i  =  1  (3.a) 

where  r,(f)  and  n,(<)  are  the  p  X  1  vectors 

*.T(0  =  l  MO  •  •  •  M<  -  [p-l)D  )]  (3-b) 


*»,r(0  - 1  MO  •  •  •  M*  -  (p-ip  )] 


s(/)  is  the  d  X  1  vector 

»r(0  =  l«i(0  • *  •  *<(01 
and  A  is  the  p  X  d  matrix 
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aile 


aue 
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(3.c) 


(3-d) 


(3.e) 


Stacking  the  p  X  1  vectors  r,(f)  (  i  =  )  into  a  mp  X  1 
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” snapshot”  vector  r{t),  we  can  further  simplify  the  notation  and  rewrite  equa¬ 
tion  (3)  as 

r(/)  =  As(i)+  b(<)  (4.a) 

where  r(t)  and  n(/)  are  the  mp  X  1  vectors 

*r(0  =  1*^(0  *'*  rm(01  (4b) 

nr(0  =  !®ir(0  '  *  *  ®m(01  (4  c) 

and  A  is  the  mp  X  d  matrix 

A0]  ane-iu'T"  •••  aut^ 

A  =  =  !  (4.d) 

Am  a  (P-»)i>  ...  a^e-jUiTm4-jvt(p-l)D 

Note  that  each  column  of  A  is  associated  with  a  different  source.  We  shall 
denote  these  column  vectors  by  A||W|  (  k  =  ),  and  refer  to  them  as  the 

direction-frequency  vectors  of  the  sources. 

Multiplying  (4. a)  by  its  conjugate  transpose  and  taking  expectations,  assum¬ 
ing  that  the  noises  are  zero-mean  with  variance  a3,  that  their  correlation  time  is 
smaller  than  D,  and  that  they  are  uncorrelated  with  each  other  and  with  the 
source  signals,  we  obtain 

R  =  ASA*  +  ft  (5. a) 

where  t  denotes  the  conjugate  transpose,  R  is  the  covariance  matrix  of  the 
received  signals,  I  is  the  identity  matrix,  and  S  is  the  covariance  matrix  of 
the  sources,  i.e., 


(5  b) 


S  =  E[  ■(!)«'(<)  |  • 

where  E  denotes  the  expectation  operator. 

We  shall  assume  that  S  is  nonsingular.  That  is,  no  source  is  a  scaled  and 
delayed  version  of  any  other  source;  in  other  words,  no  two  sources  are  fully 
correlated.  We  do,  however,  allow  for  sources  to  be  partially  correlated,  in  the 
sense  that  the  covariance  matrix  of  the  sources  may  be  nondiagonal,  as  long  as  it 
is  nonsingular.  This  allowance  makes  the  model  applicable  to  certain  kinds  of 
multipath  problems,  e.g.  those  where  the  reflection  introduces  some  random  per¬ 
turbation  to  the  multipath  signal. 

B.  The  Wideband  Problem 

Consider,  as  before,  an  array  with  m  sensors.  Assume  now  that  d  wide¬ 
band  sources  (<f  <  m),  with  identical  bandwidth  B,  impinge  on  the  array  from 
directions  9x,...,6d.  The  signals  emitted  by  the  sources  are  assumed  to  be  station¬ 
ary  zero-mean  stochastic  processes.  Using  the  notation  of  (1),  the  signal  received 
at  the  i-th  sensor  can  be  expressed  as 

n(0=  E  «*«*(*  -ra)  +  n<(0  (6) 

k-i 

Unlike  the  narrowband  problem,  where  we  have  formulated  the  problem  in 
terms  of  the  sampled  data,  it  turns  out  that  for  the  wideband  problem  it  is  more 
convenient  to  formulate  the  problem  in  terms  of  the  continuous  signals.  Assum¬ 
ing  that  we  are  observing  the  received  signals  over  a  finite  interval  T,  it  follows 
that  we  can  approximately  represent  the  received  signal  r,()  by  a  Fourier-series 


I+M 


where  are  the  Fourier-coefficients,  given  by 


*<“•'-  Ft  d‘ 


w.  =  ——  n  n  =  +  M 


•where  and  W/+jj#  are  the  lowest  and  highest  frequencies,  respectively, 

included  in  the  bandwidth  B.  Note  that  there  are  only  positive  frequencies  since 
we  are  using  the  complex  (analytic)  signal  representation,  and  that  there  are  a 
finite  number  of  them  since  we  are  assuming  that  the  signals  have  an  approxi¬ 
mately  defined  finite  bandwidth  B. 

Then  representing  both  sides  of  (6)  by  their  Fourier-coefficients,  assuming 
that  the  observation  time  is  much  longer  than  the  propagation  time  of  the  signals 
across  the  array  (r{k  «  T)  so  that  to  a  good  approximation  the  time-delay 
transforms  to  phase-shift  in  the  Fourier-domain,  we  obtain 

=  E  sk( u/„)  +  N.-K)  (8) 

*-i 


In  matrix  notation  this  becomes 


R(WJ  =  A(WJSK)+  N(W|) 


where  R(u;„)  and  N(w„)  are  the  m  X  1  vectors 
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rTK)  =  [RXM  ■  •  •  Rm{un)} 

Nr(W,)  =  I N.M  •••  NmM) 

S(u/„)  is  the  d  X  1  vector 

s  T(^«)  *  [Si(w,)  ^(w,)] 

and  A(ww)  is  the  m  X  d  matrix 
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(9.b) 


(9  c) 


(9.d) 


(9.e) 


Note  that  each  column  of  A  (a;,)  is  associated  with  a  different  source. 

We  shall  denote  these  column  vectors  by  Aft(u>#)  (  k  =  1 . d  ),  and  refer  to 

them  as  the  direction  vectors  of  the  sources. 

Multiplying  (9.a)  by  its  conjugate  transpose,  and  taking  expectations,  assum¬ 
ing  that  the  noises  are  zero-mean  and  uncorrelated  with  the  signals,  we  obtain 

£|RK)BV,)l=AK)£lSK)SlK)lAt(w.)+  £INh)N'Kl  (10) 

Next,  assuming  that  the  observation  time  is  large  compared  to  the  correlation 
time  of  the  processes  involved,  the  covariance  matrix  of  the  Fourier-coefficient 
vector  will  be  approximately  equal  to  the  spectral  density  matrix  (see,  e.g., 
Whalen  (1971)  p.81).  Thus,  assuming  that  the  noises  are  uncorrelated  with  each 
other  and  have  the  same  spectral  densities,  we  can  rewrite  (10)  as 
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Kk)  =  Ak)P(W|)Atk)  +  **( u,.)  I  (11) 

where  K(uB)  and  P{u>B)  are  the  spectral  density  matrices  of  the  processes 
{*■,(  )}, and  {«,•(  )}/_!,  respectively,  and  is  the  spectral  density  of  the 

noises 

We  shall  assume  that  P(wB)  is  nonsingular,  i.e.,  no  source  is  a  scaled  and 
delayed  version  of  any  other  source,  s  in  the  narrowband  case,  we  allow  for 
sources  to  be  partially  correlated,  in  the  sense  that  the  spectral-density  matrix  of 
the  sources  may  be  nondiagonal. 


ID.  Time-domain  Eigenstructure  Method  for  the  Narrowband  Problem 

The  solution  we  are  going  to  present  to  the  narrowband  problem  is  based  on 
the  the  eigenstructure  of  the  covariance  matrix  R  and  therefore  is  referred  to 
as  a  time-domain  method.  It  presumes  that  the  following  conditions  hold: 

(I)  The  covariance  matrix  of  the  sources  S  is  nonaingular. 

(II)  Any  set  of  d+  1  direction-frequency  vectors  are  linearly  independent. 

With  these  assumptions  it  can  be  easily  verified  that  the  eigenvalues  and 
eigenvectors  of  R,  denoted  by  {Xj  >  X2  .  •  -  >  Xmp}  and 
{V i,V2,  .  .  .  ,  V mp } ,  respectively,  have  the  following  properties: 

The  minimal  eigenvalue  of  R  is  o1  with  multiplicity  mp  -  d,  i.e.  , 

X^+ 1  =  Xrf+  2  =  •  *  •  =  Xmp  =  a2  (12) 

The  eigenvectors  corresponding  to  the  minimal  eigenvalue  are  orthogonal  to  the 
columns  of  the  matrix  A ,  i.e.  , 

A^V,-  =0  i  =  d+  l,...,mp  (13.a) 

or,  more  explicitly 

O^rf+i  >  •  •  •  »Vmp)  -L  {A||Wl ,...,  A  f<w<)  (13.b) 

where  A#Wi  is  the  »-th  column  of  A,  referred  to  as  the  direction-frequency 
vector  of  the  i-th  source. 

Note  that  for  properties  (12)  and  (13)  to  hold,  it  is  only  required  that  d, 
rather  than  d  +  1  as  in  condition  (II)  above,  direction-frequency  vectors  be 


-13- 


linearly  independent.  However,  the  stronger  condition  (II)  is  needed  to  assure  the 
uniqueness  of  the  solution  obtained  by  the  eigenstructure  method,  as  we  shall 
presently  show. 

The  eigenstructure  method  is  based  on  straightforward  exploitation  of  pro¬ 
perties  (12)  and  (13).  Observe,  first,  that  it  follows  from  (12)  that  the  number  of 
sources  can  be  determined  from  the  multiplicity  of  the  smallest  eigenvalue. 
Second,  it  follows  from  the  orthogonality  relation  (13)  between  the  direction- 
frequency  vectors  of  the  impinging  sources  and  the  eigenvectors  corresponding  to 
the  minimal  eigenvalue,  that  the  directions-of-arrival  and  the  center-frequencies  of 
the  sources  can  be  determined  simply  by  searching  for  those  direction-frequency 
vectors  that  are  orthogonal  to  the  eigenvectors  corresponding  to  the  minimal 
eigenvalue.  Note  that  to  avoid  ambiguities,  no  direction-frequency  vector  other 
than  those  corresponding  to  the  impinging  sources,  should  be  orthogonal  to  the 
eigenvectors  corresponding  to  the  minimal  eigenvalue.  That  will  be  the  case  if 
the  span  of  the  d  direction-frequency  vectors  of  the  impinging  sources  does 
not  include  any  other  direction-frequency  vector.  This  implies  that  to  assure 
unambiguous  results,  any  set  of  rf+l  direction-frequency  vectors  should  be 
linearly  independent,  as  stated  in  condition  (II)  above. 

A.  Determination  of  the  Number  of  Sources 

The  method  for  determining  the  number  of  sources  we  have  outlined  above 
was  based  on  the  eigenvalues  of  the  true  covariance  matrix.  In  practice,  however, 
the  true  covariance  matrix  is  unknown.  To  apply  the  method  we  must,  therefore, 


estimate  the  eigenvalues  from  the  data.  The  problem  is  that  the  estimated  eigen¬ 
values  will  not  obey  relation  (12).  With  probability  one,  the  small  mp-d  eigen¬ 
values  will  all  be  different,  somehow  clustered  around  their  true  value.  Thus  it  is 
impossible  to  determine  the  number  of  sources  simply  by  observing  the  small 
eigenvalues,  and  a  more  sophisticated  procedure,  based  on  statistical  considera¬ 
tions  is  needed. 

Such  a  procedure  was  developed  by  Bartlett  and  L&wley  (Bartlett  (1954), 
Lawley  (1956))  in  the  context  of  factor  analysis.  The  Bartlett-Lawley  procedure 
takes  the  form  of  a  sequence  of  nested  hypothesis  tests  for  the  multiplicity  of  the 
smallest  eigenvalue,  starting  from  the  highest  possible  multiplicity  and  testing 
successively  down  to  the  lowest  possible  multiplicity.  At  the  k’th  step  of  the  pro¬ 
cedure  (  k  ==  0,...,mp-l  ),  the  null  hypothesis  that  the  smallest  eigenvalue  of  the 
covariance  matrix  has  multiplicity  mp  -  k,  namely, 


Hk  '•  * *+ 1  =  ^*+2  '  '  ‘  = 

is  tested.  The  likelihood-ratio  for  this  problem,  under  Gaussian  assumptions,  is 
given  by 


<?*== 


ff  i, 

-*+ 1 


,<■=77  2  '.) 
[  mP~k  ,«*+  1 


N 


(14) 


where  /,  >  l2  *  *  •  >  lmp 
matrix  R,  defined  by 


are  the  eigenvalues  of  the  sample- covariance 
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R  =  4r  E  (IS) 

iV  i-l 

Note  that  Qk  is  a  monotonic  function  of  the  ratio  of  the  arithmetic- mean 
and  the  geometric-mean  of  the  smallest  m  -  k  eigenvalues  whose  equality  is 
tested. 

The  exact  distribution  of  the  statistic  Qk  under  the  null  hypothesis  Hk  is 
unknown.  However  asymptotically,  as  N— ►  oo,  it  follows  from  the  general  theory 
of  likelihood-ratios  (see  e.g  Cox  and  Hinckly  (1974))  that  the  statistic  -21nQ*  has 
a  x2  distribution  with  (mp-kf  -  1  degrees  of  freedom.  Therefore,  if  N  is 
large  an  approximate  test  of  size  a  of  Hk  is  to  reject  Hk  if 
-21n(?4>c(a;(mp-/t)2-l),  where  c(o;r)  is  the  upper  lOOo  %  point  of  the  x2 
distribution.  Note  that  the  size  a  of  the  test  is  a  parameter  left  to  the  subjec¬ 
tive  decision  of  the  designer. 

The  hypotheses  Hk  {k  =  0,l,...,mp-l)  are  tested  sequentially.  The  value 
of  k  for  which  Hk  is  first  accepted  is  selected  as  the  estimate  d  of  the  number 
of  sources.  For  the  1-D  problem  (p=l)  this  procedure  reduces  to  that  presented 
by  Simkins  (1980). 

We  should  note  that  a  different  approach  to  the  problem,  based  on  applying 
the  information  theoretic  criteria  for  statistical  model  identification  introduced  by 
Akaike,  Schwartz  and  Rissanen,  was  recently  presented  by  Wax  and  Kailath 
(1983b).  In  this  approach  no  subjective  judgement  is  required  in  the  decision  pro¬ 


cess. 
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B.  Estimation  of  the  Eigenstructure 


The  method  we  have  outlined  for  the  determination  of  the  directions-of- 
arrival  and  the  center-frequencies  was  based  on  the  eigenstructure  of  the  true 
covariance  matrix.  The  problem  is  that  in  practice  the  true  covariance  matrix 
and  hence  its  eigenstructure  are  unknown.  To  apply  the  method  we  must,  there¬ 
fore,  estimate  the  eigenstructure  from  the  data.  To  this  end,  let  the  eigenvalues 
and  eigenvectors  of  the  sample-covariance  matrix  R ,  defined  in  (15),  be  given 
by  >  /2  >...  >  lmp  and  C,,  .  .  .  ,Cmf>,  respectively.  Now,  using  the  esti¬ 
mate  of  the  number  of  sources  d,  it  follows  from  Anderson  (1963)  that  the  max¬ 
imum  likelihood  estimates  of  the  eigenvalues  and  eigenvectors  of  R  are  given 
by 


X,  =  /. 


*'  =  1,-f* 


(16.a) 


(16  b) 


V 


»  =  1  ,..,d 


(16  c) 


and  up  to  an  orthogonal  transformation  from  the  right  (the  eigenvectors 
corresponding  to  multiple  eigenvalue  define  a  subspace  and  are  nonunique  up  to 
an  orthogonal  transformation) 

V4  =  C,  is  d+ 1  ,...,mp  (16.d) 
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C.  Estimation  of  the  Directions- Of- Arrival  and  Center-Frequencies 


With  estimates  of  the  number  of  sources  and  the  eigenvectors  we  now 
proceed  to  the  simultaneous  estimation  of  the  directions-of-arrival  and  the 
center-frequencies  of  the  sources.  We  have  seen  that  the  direction-frequency  vec¬ 
tors  corresponding  to  the  d  sources  are  orthogonal  to  the  true  eigenvectors 

V,+  i . Vmp  corresponding  to  the  repeated  smallest  eigenvalue.  However, 

since  we  only  have  estimates  of  these  eigenvectors,  the  orthogonality  relation  does 
not  hold  any  more.  Instead,  the  cosine  of  the  angle  between  each  of  the 
direction-frequency  vectors  A<lWl,  .  .  . ,  A|JU<  and  each  of  the  eigenvectors 
Vj+  . ,  Vmp  will  probably  be  "close"  to  zero,  that  is 


cos  =  I  A/tW/Vj  |2  /  |A/tW|A,4_4 1  0  , 


* 

* 


l,..,d 

d+  l,...,mp 


(17) 


where  |  •  |  denotes  the  modulus  of  the  complex  number.  Thus,  if  A  t  u  denotes 

the  direction-frequency  vector  corresponding  to  a  source  at  direction-of- arrival  0 

* 

and  center-frequency  u,  then  the  d  sources  should  be  chosen  as  those  whose 
direction-frequency  vectors  are  "most  nearly  orthogonal"  to  the  set  of  eigenvec- 
tors  {Vj+  j, . . .  t  Vap }  corresponding  to  the  minimal  eigenvalue. 

There  are  many  metrics  that  can  be  applied  to  measure  this  "distance  from 
orthogonality”.  One  metric  is  the  arithmetic-mean  of  the  square  cosines  of  the 
angle  between  the  direction-frequency  vector  and  the  eigenvectors  corresponding 

to  the  smallest  eigenvalue,  that  is,  — | A/^Vj  |2/| A/„A#J . 
Another  metric  is  the  geometric-mean  of  these  quantities,  that  is 


2 

rf  iA;„vir»-'/|A/„A.,j.  Both  these  metrics  are  two  special  cases  of 
i-J+ 1 

the  general  family  (see  e.g  Pisarenko  (1072))  given  by 

i 

( — r  |ALVi|2r)r/|A;wA#w|.  Clearly,  the  arithmetic  -mean 

TOP-rf,-i+  i 

corresponds  to  r  =  1.  The  geometric-mean  can  be  shown  to  correspond  to  the 
limit  obtained  as  r-*  0. 

With  these  "orthogonality  metrics”,  we  can  now  proceed  to  the  estimation  of 
the  direction-of-arrival  and  the  center-frequencies  of  the  sources.  Note  that  since 
we  are  interested  only  in  the  extremal  points  of  these  metrics,  we  can  extract  the 
information  from  any  monotonic  function  of  them.  In  order  to  resemble  the  form 
of  Capon’s  Minimum  Variance  (MV)  and  Burg’s  Maximum  Entropy  (ME)  estima¬ 
tors,  it  seems  natural  to  choose  the  inverse  of  these  metrics  as  the  representative 
form.  For  the  arithmetic-mean  metric,  we  therefore  plot 

m ->  =  ,  -  <»»..) 

— *"7  2  ia/.v.i* 

For  the  geometric-mean  metric,  we  plot 

m  W)  = - •A;"'A,“I  l  (18  b) 

ff  |A/„V,|-»-* 

I-J+ 1 


and  for  the  general-mean  metric,  we  plot 
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MO*  U)  =  - 7  (18.C) 

t — *-j  £  ialv,ij')' 

The  directions-of-arrival  and  the  center-frequencies  of  the  impinging  sources 
are  determined  as  the  c!  points  (0i,u?i)  (0j,wj)  that  yield  the  highest 

peaks  of  either  of  these  2-D  functions. 

The  estimators  presented  above  differ  in  their  resolution  and  accuracy  pro¬ 
perties.  The  arithmetic-mean  estimator  will  have  a  peak  in  the  point  {0,  u)  if 
and  only  if  A  is  "almost  orthogonal”  to  all  the  eigenvectors 
The  geometric-mean  estimator,  on  the  other  hand,  will  have  a 
significant  peak  at  the  point  (9,  u>)  even  if  Atu  is  "almost  orthogonal”  to  only 
one  of  the  eigenvectors  V^+I>  .  . .  jV^p.  Therefore,  the  geometric-mean  esti¬ 
mator  will  have  higher  resolution  but  lower  accuracy  than  the  arithmetic- mean 
estimator. 

More  delicate  trade-offs  between  resolution  and  accuracy  can  be  obtained  by 
using  the  general  estimator  (18.c)  with  the  parameter  r  chosen  appropriately. 

We  should  note  that  Sharman  et.  al  (1683)  have  recently  shown  that  the 
arithmetic-mean  metric  (18. a)  is  asymptotically  optimal  for  the  case  of  a  single 
source.  That  is,  under  a  small  error  assumption,  the  bias  of  the  estimator  (18.a) 
approaches  zero,  and  its  variance  approaches  the  Cramer-Rao  lower  bound. 

D.  Estimation  of  the  Covariance  Matrix  of  the  Sources 


The  source  covariance  matrix  S  contains  valuable  information  on  the 


impinging  sources.  The  diagonal  elements  indicate  the  power  of  the  sources,  and 
the  off-diagonal  elements  indicate  the  cross-correlation  between  them.  This  infor¬ 
mation  can  be  valuable  for  distinguishing  between  a  direct  path  and  multipath, 
for  example. 

Using  the  estimates  of  the  number,  the  directions-of-arrival  and  the  center- 
frequencies  of  the  sources,  the  estimate  of  the  covariance  matrix  of  the  sources 
follows  immediately  from  the  works  of  Schmidt  (1079)  and  of  Bienvenu  and  Kopp 
(1981).  Note  first  that  (12)  implies  that 

ASA'-  £  (X.  -^V.Vt  (10) 

1—1 

Solving  this  equation  for  S,  and  substituting  estimated  quantities  for  the  true 
quantities,  yields 


where 


S  =  (A  *A  )_,A  *V  AV^A  (A  *A  )-1 


V 


^ 

X,  -  d2  0 

0  X*  -  d* 


(20.a) 


(20.b) 


(20.c) 


(20.d) 


with  denoting  the  direction-frequency  vector  corresponding  to  the 


estimated  direction-of-arrival  and  center- frequency  of  the  i’th  source. 
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E.  Discussion 

The  maximum  number  of  sources  that  the  method  can  handle  is  pm  -  1, 
assuming  that  not  more  than  m  -  1  of  them  are  “co-frequency"  (i.e.,  have  the 
same  center-frequency).  These  restrictions  stem  from  condition  (D)  of  our  under¬ 
lying  assumptions  which  requires  that  any  set  of  d  +  1  direction-frequency 
vectors  should  be  linearly  independent.  Thus,  since  the  frequency-direction  vector 
is  of  dimension  mp,  the  number  of  sources  d  must  be  no  greater  then  pm  -  1. 
The  restriction  on  the  number  of  co-frequency  sources  is  because  the  subspace  of 
all  possible  co-frequency  direction-frequency  vectors  is  of  dimension  no  greater 
than  the  the  number  of  sensors  m,  regardless  of  the  number  of  taps  p.  That  is, 
any  set  of  m  or  more  co-frequency  direction-frequency  vectors  will  definitely 
be  linearly  dependent.  Therefore,  for  condition  (II)  above  to  hold,  the  number  of 
co-frequency  sources  should  be  no  greater  than  m  -  1. 

The  structure  of  the  array  for  which  the  method  is  applicable  is  also  res¬ 
tricted  by  condition  (II)  above.  This  condition  does  not  hold  in  general  for  arbi¬ 
trary  geometry  and  arbitrary  number  of  taps.  It  turns  out,  for  example,  that  for 
the  case  of  a  uniform  linear  array,  this  condition  does  not  hold  for  one  tap,  i.e., 
p  =  1.  This  is  because  all  the  direction-frequency  vectors  corresponding  to 
points  (0  u/)  on  the  curve  w  sin#  =  constant  are  identical,  and  hence 
linearly  dependent.  Thus  to  assure  the  applicability  of  the  method  in  the  case  of 
a  uniform  linear  array,  a  tapped-delay-line  with  more  than  one  tap  is  needed. 
The  exact  number  of  taps  p  that  will  guarantee  the  applicability  of  the  method 
for  a  given  geometry  and  frequency  band  is  unknown  even  for  the  simple  case  of 
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uniform  linear  array.  Nevertheless,  the  simulation  results  (see  sec.  V)  show  that 
this  number  is  usually  very  small. 

The  method  presented  (with  the  arithmetic-mean  metric)  reduces  to 
Pisarenko  and  Schmidt’s  method  for  the  1-D  problem,  that  is,  for  the  case  that 
all  the  sources  have  the  same  known  center-frequency  and  only  one  tap  is  used 
(p  =  1)  (Pisarenko's  method  is  a  special  case  of  Schmidt’s  method  where  only 
one  of  the  eigenvectors  corresponding  to  the  smallest  eigenvalue  is  used). 
Larimore’s  (1081)  method  and  the  variant  of  Reddy  et  al.  (1082b)  essentially 
amounts  to  on-line  computation  of  the  eigenvector  corresponding  to  the  smallest 
eigenvalue  of  the  covariance  matrix  (a  fact  not  pointed  out  in  either  paper). 
Thus,  Their  method  can  be  regarded  as  an  on-line  implementation  of  the  off-line 
2-D  extension  to  Pisarenko’s  method  we  have  presented  above. 


IV.  Frequency- domain  Eigenstructure  Method  for  the  Wideband  Problem 

The  solution  we  are  going  to  present  for  the  wideband  problem  is  based  on 
the  the  eigenstructure  of  the  spectral-density  matrix  K(w„)  and  therefore  is 
referred  to  as  a  frequency-domain  method.  It  presumes  that  the  following  condi¬ 
tions  hold: 

(I)  The  spectral-density  matrix  P(u;)l)  is  nonsingular. 

(II)  Any  set  of  d+  1  direction  vectors  are  linearly  independent. 

With  these  assumptions  it  can  be  easily  verified  that  the  eigenvalues  and 
eigenvectors  of  K  (<*>„),  denoted  by  Xjfa/,,)  <  •  •  <  Am  (u^„ )  and 
V  j(u;„  m(w„)  ,  respectively,  have  the  following  properties: 

The  minimal  eigenvalue  of  K(u;,)  is  <r2(u >„)  with  multiplicity  m  -  d,  i.e.  , 

xrf+lK)  =  •••  =  Xm(wJ  =  (21) 

The  eigenvectors  corresponding  to  the  minimal  eigenvalue  are  orthogonal  to  the 
columns  of  the  matrix  A(w.),  i.e.  , 

Atwjty.fwj  =  0  i  =  d+ 1 (22. a) 

or,  more  explicitly 

(v  i(w.),...,V  m(w„)}  -L  {Ati(uin)  A  #<  (<*;„)}  (22. b) 

where  A|((wJ  is  the  it  ’th  column  of  the  matrix  A  (a;,),  referred  to  as  the 
direction-vector  of  the  k  ’th  source.  Note  that  the  orthogonality  relation  (22) 
holds  for  every  un  that  is  included  in  the  signal  bandwidth  B. 
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These  relations  form  the  basis  for  the  frequency-domain  eigenstructure 
method.  Because  of  the  analogy  with  the  time-domain  method  we  shall  proceed 
directly  to  the  description  of  the  method. 


A.  Determination  of  the  Number  of  Sources 


As  iu  the  time-domain  method,  the  number  of  sources  can  be  determined 
from  the  multiplicity  of  the  smallest  eigenvalue  of  the  spectral  density  matrix 
K(w„).  Applying  the  Bartlett-Lawley  procedure,  the  hypothesis  tested  at  the  k 
’th  step  is  given  by 


^i(wn)  •  ^*+l(Wn)  ...  — 


Following  a  well-known  analogy  between  multivariate  analysis  and  time-series 
analysis  (see,  e.g,  Brillinger  (1964)  and  Wahba  (1968)),  the  likelihood-ratio  for 
this  hypothesis  is  obtained  simply  by  substituting  the  eigenvalues  of  the  periodo- 
gram  estimate  of  the  spectral-density  matrix  for  those  of  the  sample- covariance 
matrix,  that  is 


QkM  = 


.  n  *(«.) 

i—k+  1 


(-TT  £ 
m~k  •-*+! 


(23) 


where  /i(wn)>/2(w»)>  '  *  *  >/m(w #)  denote  the  eigenvalues  of  the  periodogram 
estimate  of  the  spectral  density  matrix 


£  R,K)R,'K) 

L  i-I 


(24) 


with  11,(0*,,)  denoting  the  Fourier-coefficient  vector  of  the  »-th  subinterval  at 
the  frequency  un  and  L  denoting  the  number  of  subintervals. 

Since  the  sources  are  wideband,  the  number  of  sources  can  be  estimated 
from  any  frequency  w,  £  B.  Thus,  to  minimize  the  probability  of  error,  the  fol¬ 
lowing  composite  hypothesis,  which  captures  the  information  in  all  the  frequen¬ 
cies  £  B,  should  be  tested: 

Hk  ■  x*+  iK)  =  .  .  •  =  Xm(wn)  for  everyo;,  6  B 

Assuming  that  the  observation  time  T  is  much  larger  than  the  processes  corre¬ 
lation  time,  it  follows  (see,  e.g,  Whalen  p.81)  that  the  Fourier-coefficients 
corresponding  to  different  frequencies  are  uncorrelated.  Thus,  under  Gaussian 
assumptions,  the  composite  likelihood-ratio  test  is  given  by 

Qk  ='nW<?*(u;n)  (25) 

The  statistic  -2 In  Qk  has  in  this  case  an  asymptotic  x2  distribution  with 
M[(m-k )2  -  1]  degrees  of  freedom.  Therefore,  if  N  is  large  an  approximate  test 
of  size  a  of  Hk  is  to  reject  Hk  if  -21nQ*>c(a;(m-fc)2-l),  where  e(a;r)  is  the 
upper  100a  %  point  of  the  \r  distribution. 

The  hypotheses  Hk  ( k  =  0,l,...,m-l)  are  tested  sequentially.  The  value  of 

* 

k  for  which  Hk  is  first  accepted  is  selected  as  the  estimate  d  of  the  number  of 
sources.  For  the  1-D  problem  (M=l)  this  procedure  reduces  to  that  presented  by 
Ligget  (1972). 

We  should  note  that  a  different  approach  to  the  problem,  based  on  applying 
the  information  theoretic  criteria  for  statistical  model  identification  introduced  by 
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B.  Estimation  of  the  Eigenstructure 

Let  ‘  ‘  •  >LM  “<*  C,(w,)(  .  .  .  ,Cm(<Jn)  denote  the 

eigenvalues  and  eigenvectors,  respectively,  of  the  the  spectral  density  matrix 
K(ua).  Following  again  the  duality  between  multivariate  analysis  and  time- 
series  analysis,  it  follows  that  the  maximum  likelihood  estimates  of  the  eigenvalue 


and  eigenvectors  of  K  (w, )  are  given  by 

M.  )«*<".)  <  =  J i  <2M 

•W--*-?  £  /,K)  (26.b) 

m~d  ,-J+| 

V,K)-C,(«.)  «=l,...,rf  (26.c) 

and  up  to  an  orthogonal  transformation  from  the  right 

,-  =  d+l . m  (26. d) 


C.  Estimation  of  the  Directlons-Of- Arrival 

By  analogy  with  the  time  domain  method,  the  d  directions-of-arrival 
should  be  chosen  as  those  whose  direction  vectors  are  "most  nearly  orthogonal” 
to  the  set  of  eigenvectors  (Vk(o;B),  k  =  l,...,m  -  d  }.  Note  that  this  "distance 
from  orthogonality”  should  be  measured  for  all  the  frequency  bins  wa  6  B. 
Thus,  we  must  first  measure  the  "distance  from  orthogonality”  at  each  frequency 
bin,  and  then  combine  the  resulted  measures  for  the  different  frequencies.  Using 
the  arithmetic-mean  metric  both  for  the  individual  frequency  bins  and  for  the 
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combination  over  the  frequency  range,  yields  the  estimator  (the  subscript  11 
stands  for  the  double  use  of  the  arithmetic-mean  metric),  given  by 

Jim  = 

Using  the  arithmetic-mean  metric  for  the  individual  frequency  bins  and  the 
geometric-mean  metric  for  the  combination  over  the  frequency  range,  yields  the 
estimator 


_ |aJ(u.)a,k)| _ 

]  1+1/  |  m  . 

S  |A/K)vkK)|* 


(27. a) 


j  io(*)  = 


_ |AjK)A,(aQl _ 

im-H  t  iA;K)Vk(o,.)i2)“ 


(27  b) 
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The  directions-of-arrival  of  the  sources  are  then  determined  as  the  d 
points  0j,  .  .  .  ,  0}  that  yield  the  highest  peaks  of  either  one  of  the  estimators 
given  in  (27). 

Note  that  other  combinations  of  the  arithmetic-mean  and  the  geometric- 
mean  metrics,  as  well  as  other  combinations  of  different  metrics  from  the  general 
family  of  metrics  introduced  in  section  ID  can  be  used.  From  the  corresponding 
discussion  for  the  time-domain  method  it  follows  that  the  estimator  (27.b)  exhi¬ 
bits  higher  resolution  and  lower  accuracy  than  the  estimator  (27.a). 


D.  Estimation  of  the  Spectral  Density  Matrix  of  the  Sources 

The  spectral  density  matrix  of  the  sources  P(w„)  contains  valuable  infor¬ 
mation  on  the  impinging  sources.  The  diagonal  elements  yield  the  spectral 
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densities  of  the  sources,  and  thus  provides  a  too!  for  classifying  the  sources 
according  to  their  spectral  ’’signatures”.  The  off-diagonal  elements  indicate  the 
amount  of  correlation  existing  between  the  sources,  and  thus  provide  a  way  to 
distinguish  between  a  direct  path  and  multipath,  for  example. 

The  estimation  scheme  is  analogous  to  that  described  for  the  time-domain 
method.  Therefore  by  analogy  with  (20),  we  obtain 

PK)=(A'K  )A(w.  )]-'A'K  )VK  |i|U.|vV.  )A(u„  )|A'(u.  )AK  )|-'(28.a) 
where 

Vk)  =  tV,K)  •••  VjK)]  (28.b) 

AK)=  ‘  .  (28. c) 

and 

A(w„)  =  [A/((wJ  •••  A#<(wJ)  (26.d) 

with  Aj(<jJn)  denoting  the  direction-vector  corresponding  to  the  estimated 
direction-of-arrival  of  the  i-th  source. 

E.  Discussion 

The  maximum  number  of  sources  the  method  can  handle  is  m  -  1.  From 
the  requirement  that  any  set  of  d  +  1  direction-vectors  be  linearly  indepen- 


dent,  it  follows  that  since  the  direct  ion-vector  is  of  dimension  m,  the  number  of 
sources  should  be  no  greater  than  m  -  1. 

The  structure  of  the  array  for  which  the  frequency-domain  method  is  appli¬ 
cable  is  restricted  by  condition  that  any  d  +  1  direction-vectors  be  linearly 
independent.  Note,  however,  that  for  a  omnidirectional  (ati  =  1)  uniform 
linear  array,  this  condition  is  always  satisfied,  because  of  the  Vandermonde  type 
of  the  direction-vectors. 

The  estimation  of  the  spectral-density  matrix  of  the  received  signals  is  by  no 
means  restricted  to  the  peridogram  method.  Any  multivariate  spectral  estimation 
technique,  parametric  or  nonparametric  can  be  used  if  there  is  reason  to  believe 
that  it  will  give  bettter  results. 

The  frequency-domain  eigenstructure  method  we  have  presented  reduces  to 
Bienvenu  and  Kopp’s  method  for  the  1-D  problem,  that  is,  when  the  sources 
occupy  only  one  frequency  bin  (M=l). 
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V.  Simulation*  Results 

In  this  section  we  present  computer  simulation  results  that  illustrate  the  per¬ 
formance  of  the  proposed  methods.  All  the  examples  will  refer  to  a  uniform  linear 
array.  For  this  type  of  array  it  is  convenient  to  use  the  notion  of  normalized  fre¬ 
quency,  defined  as  -  (  where  w  is  the  frequency  of  the  source,  0  is  its 

direction-of-arrival,  A  is  the  spacing  between  the  sensors,  and  c  is  the  speed  of 
propagation. 

In  the  first  example  we  wanted  to  demonstrate  the  time-domain  method  for 
the  narrowband  problem.  We  considered  2  sinusoidal  sources,  having  normal¬ 
ized  frequencies  0.2  and  0.3  and  normalized  wavenumbers  0.125  and  0.2, 
respectively,  impinging  on  a  uniform  linear  array  of  3  sensors  (m=3),  each  fol¬ 
lowed  by  a  tapped-delay-line  with  3  taps  (p=3).  The  signal-to-noise  ratio  was 
lOdB.  The  results,  obtained  from  200  "snapshots”,  using  the  estimator  (18.a),  are 
presented  in  Fig.  2.  The  two  peaks  corresponding  to  the  two  sources  are  clearly 
seen. 

In  the  second  example  we  wanted  to  demonstrate  the  ability  of  the  time- 
domain  method  to  resolve  more  sources  than  the  number  of  sensors,  given  that 
the  number  of  taps  is  appropriate.  We  considered  4  sinusoidal  sources,  having 
normalized  frequencies  0.1  ,  0.2  ,  0.3  and  0.4  and  normalized  wavenumbers 
0.125 , 0.25 , 0.35  ,  and  0.42  ,  respectively,  impinging  on  a  uniform  linear  array 
of  3  sensors  (m=3),  each  followed  by  a  tapped-delay-line  of  3  taps  (p=3).  The 
signal-to-noise  ratio  was  lOdB.  The  results,  obtained  from  1000  "snapshots", 


using  the  estimator  (18.a),  are  presented  in  Fig.3.  The  four  peaks  corresponding 
to  the  four  sources  are  clearly  seen. 

In  the  third  example  we  wanted  to  demonstrate  the  frequency-domain 
method  for  the  wideband  problem.  We  considered  3  wideband  sources  having 
identical  spectra,  centered  at  0.25  with  bandwidth  0.05,  impinging  from  36°, 
60°  ,  and  120°  on  a  uniform  linear  array  of  5  sensors  (m=5)  spaced  a  third 
wavelength  apart.  The  $ignal-to-noise  ratio  was  6dB.  The  results  obtained  from 
100x64  samples,  using  the  estimator  (27.a)  are  presented  in  Fig.  4  (a  64  point 
DFT  was  used  to  compute  the  spectral-density  matrix).  The  three  peaks 
corresponding  to  the  three  sources  are  clearly  seen. 


VI.  Concluding  Remarks 


The  eigenstructure  methods  devloped  by  Pisarenko  (1073),  Schmidt  (1070) 
and  Bienvenu  and  Kopp  (1080),  were  confined  to  the  1-D  problem.  In  this  paper 
we  have  generalized  these  methods  to  the  2-D  problem,  that  is,  to  the  simultane¬ 
ous  estimation  of  the  spatio-temporal  spectrum  of  the  signals  received  by  the  pas¬ 
sive  array.  We  have  presented  both  a  time-domain  method,  based  on  the  eigen¬ 
structure  of  the  covariance  matrix  of  the  received  signals,  and  a  frequency- 
domain  method,  based  on  the  eigenstructure  of  the  spectral-density  matrix  of  the 
received  signals.  Though  the  time-domain  method  was  applied  to  the  nar¬ 
rowband  problem  and  the  frequency-domain  method  was  applied  to  the  wideband 
problem,  the  methods  are  in  fact  applicable  to  both  problems.  The  applicability 
of  the  frequency-domain  method  to  the  narrowband  case  can  be  seen  by  noting 
that  the  wideband  problem  includes  the  narrowband  problem  as  a  special  case 
corresponding  to  M  —  1.  The  applicability  of  the  time-domain  method  to  the 
wideband  problem  can  be  seen  by  noting  that  from  the  Fourier-representation  it 
follow  that  a  wideband  signal  can  be  decomposed  to  a  weighted  sum  of  A/  com¬ 
plex  exponentials.  Thus,  the  applicability  of  the  time-domain  method  to  the 
wideband  problem  is  guaranteed  if  the  number  of  taps  p  is  at  least  equal  to  M. 
However,  the  computational  complexity  involved  in  implementing  the  time- 
domain  method  for  the  wideband  problem  and  the  frequency-domain  method  for 
the  narrowband  problem,  make  these  alternatives  unattractive  in  practice. 

The  resolution  offered  by  the  eigentructure  methods  depends  crucially  on  the 
quality  of  the  estimates  of  the  covariance  and  spectral-density  matrices.  If  the 


data  record  is  long  enough  to  enable  resonably  good  estimates  of  these  matrices, 
the  resolution  is  very  high  even  in  relatively  low  signal-to-noise  ratios.  However, 
if  the  data  record  is  too  short,  the  performance  of  these  methods  deterioate  drast¬ 
ically. 

Theoretical  analysis  of  the  performance  of  these  method  for  finite  record 
length  seems  to  be  a  very  difficult  problem.  The  only  analytical  results  available 
are  for  the  asymptotic  where  the  record  length  approaches  infinity.  It  has  been 
recently  shown  by  Sharman  et  al  (1083)  that  the  eigenstructure  estimator  based 
on  the  arithmetic-mean  metric  is  asymptotically  unbiased,  and  that  for  the  case 
of  a  single  source,  it  is  even  efficient,  that  is,  its  variance  approches  the  Cramer- 
Rao  lower  bound. 

The  methods  we  have  presented  can  be  straightforwardly  extended  to  the 
case  that  the  sources  have  diverse  polarization  (see  Schmidt  (1970)  and  Ferrara 
and  Parks  (1983))  and  to  the  more  general  problem  of  source  localization  (see 
Wax  et.  al  (1982)). 

Acknowledgement 

The  authors  wish  to  acknowledge  helpful  discussions  with  Dr.  B.  Medoff. 


REFERENCES 


[1]  Anderson,  T.  W.  (1903):  "Asymptotic  Theory  for  Principal  Component 
Analysis,”  Ann.  of  Math.  Stat.,  Vol.  34,  No.  1,  1903,  pp.  122*148. 

[2]  Bienvenu,  G.  and  Kopp,  L.  (1980):  "Adaptivity  to  Background  Noise  Spa¬ 
tial  Coherence  for  High  Resolution  Passive  Methods”,  in  Proe.  IEEE 
ICA  qSP  80,  (Denver,  CO),  pp.  307*310. 

[3]  Bienvenu,  G.  and  Kopp,  L.  (1981):  "Source  Power  Estimation  Method 
Associated  with  High  Resolution  Bearing  Estimation”,  in  Proe.  IEEE 
ICASSP,  (Atlanta,  GA),  pp.  302-305. 

|4]  Brillingcr,  D.R.  (1964):  "A  Frequency  Approach  to  the  Techniques  of  Prin¬ 
ciple  Components,  Factor  Analysis  and  Canonical  Variates  in  the  Case  of 
Stationary  Series,”  Royal  Statistical  Society  Conference,  (Cardiff,  Wales). 

[5]  Bronez,  T.P.  and  Cad  row,  J.A.  (1983):  "An  Algebraic  Approach  to  super¬ 
resolution  Array  Processing,"  IEEE  Trans,  on  AES,  Vol  19,  No  1,  pp.123- 
133. 

[6]  Cantoni,  A.  and  Godara,  L.C.  (1980):  "Resolving  the  Directions  of  the 
Sources  in  a  Correlated  Field  Incidenting  on  an  Array,"  JASA,  Vol.  67, 
No.  4,  pp.  1247-1255. 

|7]  Cox,  D.R.,  and  Hinkley,  D.V.  (1974):  Theoretical  Statistics.  Chapman  and 
Hall,  London. 

|8]  Ferrara,  E.R  and  Parks,  T.M.  (1983):  "Direction  Finding  with  an  Array  of 
Antennas  Having  Diverse  Polarizations,"  IEEE  Trans,  on  AP,  Vol  31,  No 
2,  pp.231-236. 

(9]  Frost,  O.L.  and  Sullivan,  T.M.  (1979):  "High  Resolution  Two-Dimensional 
Spectral  Analysis,"  in  Proe.  ICASSP  79,  (Washington,  DC),  pp.  673-676. 

(10]  Hahn,  W.R.  and  Tretter,  S.A.  (1973):  "Optimum  Processing  for  Delay- 
Vector  Estimation  in  Passive  Signal  Array,”  IEEE  Trans,  on  IT,  Vol.  19, 
No.  5,  pp.  608-614. 

(11]  Jackson,  L.B.  and  Chien,  C.H.  (1979):  "Frequency  and  Beering  Estimation 
by  Two-Dimensional  Linear  Prediction,”  in  Proe.  IEEE  ICASSP  79, 
(Washington,  DC),  pp.  665-668. 

(12]  Jain,  A.K.  and  Raganath,  S.  (1978):  "Two-Dimensional  Spectral  Estima¬ 
tion,"  in  Proe.  RDAC  Spectrum  Estimation  Workshop,  (Rome,  N.Y),  pp. 
151-157. 


-36- 


[13]  Johnson,  D.H.  and  Degraff,  S.R.  (1982):  "Improving  the  Resolution  of 
Bearing  in  Passive  Sonar  Arrays  by  Eigenvalue  Analysis,”  IEEE  Tram,  on 
ASSP,  Vol  30,  No  4,  pp.  638-647. 

[14]  Kumaresan,  R.  and  Tufts,  D.W.  (1981):  "A  Two-Dimensional  Technique 
for  Frequency-Wavenumber  Estimation,”  Proe.  IEEE  (Lett.),  Vol  69,  No 
11,  pp.  1515-1517. 

[15]  Kumaresan,  R.  and  Tufts,  D.W.  (1983):  "Estimation  of  Arrival  of  Multiple 
Plane  Waves,"  IEEE  Tram,  on  AES,  Vcl  19,  No  1,  pp.  134-139. 

[16]  Larimore,  M.G.  (1981):  "Adaptive  Two-dimensional  Harmonic  Retrieval,” 
Proe.  15th  Atilomar  Conf.  Cir.  Sytt.  Comp.,  November  1981. 

[17]  Ligget,  W.S.  (1972):  "Passive  Sonar:  Fitting  Models  to  Multiple  Time 
Series,”  in  Nato  ASI  on  Signal  Procuring,  (Loughborough,  U.K.),  pp.  327- 
345. 

[18]  Lim,  J.S.  and  Malik,  N.L.  (1981):  "A  New  Algorithm  for  Two-Dimensional 
Maximum  Entropy  Power  Spectrum  Estimation,”  IEEE  Tram,  on  ASSP, 
Vol  29,  No  3,  pp.  401-413. 

[19]  Lawley,  D.N.,  (1956):  "Tests  of  Significance  of  the  Latent  Roots  of  the 
Covariance  and  Correlation  Matrices,”  Biometriea,  43,  pp.  128-136. 

[20]  Good,  I.J.  (1963):  "Weighted  Covariance  for  Estimating  the  Direction  of  a 
Gaussian  Source,”  in  Time  Series  Analysis,  M.  Rosenblat  Ed.,  New  York, 
Wiley,  1903. 

[21]  McClellan,  J.H.  (1982):  "Multidimensional  Spectral  Estimation,”  Proe. 
IEEE,  Vol  70,  No  9,  pp.  1029-1039. 

[22]  Morf,  M.  et  al.  (1981):  "Investigation  of  New  Algorithms  for  Locating  and 
Identifying  Spatially  Distributed  Sources  and  Receivers,”  DARPA  Tech. 
Rep.  No.  M355-1. 

[23]  Nehorai,  A.  and  Morf,  M.  (1982):  "Estimation  of  Time  Difference  of 
Arrival  for  Multiple  ARMA  Sources  by  Pole  Decomposition,"  IEEE  Trans 
on  ASSP,  Vol  31,  No  6. 

[24]  Owsley,  N.L.  (1973):  "A  Recent  Trend  in  Adaptive  Spatial  Processing  for 
Sensors  Arrays:  Constrained  Adaptation,”  in  Signal  Procuring,  (J.W.R. 
Griffiths  et  al.  Eds.),  New  York:  Academic  Press. 

[25]  Pisarenko,  V.F.  (1972):  "On  the  Estimation  of  Spectra  by  Means  of  Non¬ 
linear  Functions  of  the  Covariance  Matrix,"  Geophyt.  J.  R.  artr.  Soe.,  Vol 
28,  pp.  511-531. 


-37- 


[26]  Pisarenko,  V.F.  (1973):  "The  Retrieval  of  Harmonics  from  Covariance 
Functions,”  Geoph.  J.  Royal  Astronomical  Soc.,  1972,  pp.  347-366 

[27]  Porat,  B.  and  Friediander,  B.  (1983):  "Estimation  of  Spatial  and  Spectral 
Parameters  of  Multiple  Sources,"  IEEE  Trans  on  IT,  Vol  29,  No  3,  pp. 
412-425. 

[28]  Reddi,  S.S.  (1979):  "Multiple  Source  Location  •  A  digital  Approach,”  IEEE 
Trans,  on  AES,  Vol  15,  No  1,  pp.  95-105. 

[29]  Reddy,  V.U.,  Egardt,  B.  and  Kailath,  T.  (1982):  "Least-squares  Type 
Algorithm  for  Adaptive  Implementation  of  Pisarenko's  Harmonic  Retrieval 
Method,”  IEEE  Trans,  on  ASSP,  Vol.  30,  No.  3,  June  1982,  pp.  399-405. 

[30]  Reddy,  V.U.,  Citron,  T.K  and  Kailath,  T.  (1982):  "Adaptive  Algorithms 
for  Two-Dimensional  Harmonic  Retrieval,"  in  Proc.  16th  Asilomar  Conf. 
Cir.,  Syst.  Comp.,  (Pacific  Groove,  CA),  pp.  177-181. 

[31]  Roucos,  S.  and  Childers,  D.G.  (1980):  "A  Two-Dimensional  Maximum 
Entropy  Spectral  Estimator,"  IEEE  Trans,  on  IT,  Vol  26,  No  5,  pp.  554- 
560. 

[32]  Sharman,  K.,  Wax,  M.  and  Kailath,  T.  (1983):  "Asymptotic  Statistics  of 
the  Eigenstructure  Spectral  Analysis  Method,”  submitted  to  IEEE  Trans, 
on  ASSP 

[33]  Schmidt,  R.O.  (1979):  "Multiple  Emitter  Location  and  Signal  Parameter 
Estimation,"  in  Proc.  RADC  Spectrum  Estimation  Workshop,  (Rome, 
N.Y),  pp.  243-258. 

[34]  Schmidt,  R.O.  (1981):  "A  signal  subspace  approach  to  multiple  emitter 
location  and  spectral  estimation,”  Ph.D.  dissertation,  Stanford  University, 
November  1981. 

[35]  Schweppe,  F.C.  (1968):  "Sensor  Array  Data  Processing  for  Multiple  Signal 
Sources,"  IEEE  Trans,  on  IT,  Vol  14,  pp.  294-305. 

[36]  Simkins,  D.N.  (1980):  "Multichannel  Angle-Of-Arrival  Estimation,"  Ph.D 
Dissertation,  Stanford  University. 

[37]  Thompson,  P.A.  (1979):  "An  Adaptive  Spectral  Analysis  Technique  for 
Unbiased  Frequency  Estimation  in  the  Presence  of  White  Noise,"  in  Proc. 
13th  Asilomar  Conf.  Cir.,  Syst.  Comp.,  (Pacific  Groove,  CA),  pp.  529-533 

[38]  Wahba,  G.,  (1968):  "On  the  Distribution  of  Some  Statistics  Useful  in  the 
Analysis  of  Jointly  Stationary  Time- Series,”  Ann.  Math.  St  at.,  Vol  39,  pp. 
1849-1802. 


-38- 


[39]  Wax,  M.,  Shan,  T-J.  and  Kailath,  T.  (1982):  "Location  and  Spectral  Den¬ 
sity  Estimation  of  Multiple  Sources,”  Proe.  16th  Asilomar  Conf.  Cir.,  Syst. 
Comp.,  (Pacific  Grove,  CA),  pp.  322-326. 

[40]  Wax,  M.  and  Kailath,  T.  (1983a):  "Optimum  Localization  of  Multiple 
Sources  by  Passive  Arrays,"  IEEE  Trans,  on  ASSP.  Vol  31,  No  5,  pp. 
1210-1218. 

[41]  Wax,  M.  and  Kailath,  T.  (1983b):  "Determining  the  Number  of  Signals  by 
Information  Theoretic  Criteria,"  ASSP  Spectrum  Estimation  Workshop  II, 
(tampa,  FL),  pp.  192-196. 

[42]  Whalen,  A.  D.  (1971):  Detection  oj  Signal s  in  Noise,  Academic  Press,  New 
York. 


Figure  3.  4  sinusoidal  sources  having  normalized  frequencies  0.1,  0.2,  0.3  and  0.4  and  normalized 
wavenumbers  0.125.  0.25,  0.35  and  0.42,  respectively,  impinge  on  uniform  linear  array 
of  3  sensors  (m*3)  followed  by  a  tapped-del ay-line  of  3  taps  (P=3).  The  signal -to-nois 
ratio  was  lOdB.  The  results  were  obtained  from  1000  "snapshots",  using  estimation  (17. 


Figure  4.  3  wideband  sources  having  identical  spectra  centered  at  0.25  with  bandwidth  0.05 
from  36°,  60°,  and  120°  on  a  uniform  linear  array  of  5  sensors  (m=5)  spaced  a  th 
wavelength  apart.  The  signal -to-noise  ratio  was  6dB.  The  results  were  obtained 
64x100  samples,  using  estimator  (25, a).  (A  64  OFT  was  used  to  compute  the  estim 
spectral -densi ty  matrix. ) 
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